from pylab import *
from scipy.signal import butter, lfilter

# Calculate Butterworth filter coefficients.
# This will be used to low-pass filter a random signal
[B,A] = butter(8, 0.3)
B.shape = (B.size,1)
A.shape = (A.size,1)

# Create arrray of random numbers, use seed to always
# generate the same sequence, so that we may test.
seed(3)
x = randn(4000,1) * 800.
# Initiate state of the filter.
nfilt = max(len(A),len(B))
#~ zi = zeros((max(len(A),len(B))-1,1))
za2 = concatenate( (eye(nfilt-2), zeros((1,nfilt-2),dtype=float64)))
za = [eye(nfilt-1) - [-A[1:nfilt-1].transpose() za2]
zb = [B[1:nfilt-1].transpose() - A[1:nfilt-1].transpose() * b[0]]
zi = za / zb
# Filter in the forward direction
x1,zf = lfilter(B, A, x, axis = 0, zi = zi)
#~ x2 = flipud(x1).copy()
#~ x3,zf = lfilter(B,A,x2, axis = 0, zi = zf)
#~ x4 = flipud(x3).copy()
